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1.  BACKGROUND 


The  proposed  work  addresses  some  research  issues  arising  from  Computational  Fluid  Dynamics 
(CFD)  simulations  performed  in  support  of  the  KTA2-12  program.  This  program  is  designed  to 
evaluate  computational  technology  for  predicting  highly  separated  flowfields  for  missile  configura¬ 
tions.  Extensive  experimental  data  at  various  axial  locations  is  available  for  the  test  problem  in  the 
form  of  surface  pressure  coefficients,  and  flowfield  pitot  pressures  for  off  surface  structures.  Two 
flow  solvers,  four  turbulence  models,  and  seven  grids  were  used  in  various  combinations  to  study 
the  dependence  of  the  predictions  on  a  given  component.  One  flow  solver  consistently  yielded  pre¬ 
dictions  closest  to  the  data.  For  this  flow  solver  the  results  were  most  sensitive  to  the  grid  used.  It 
is  desirable  that  the  point  distribution  reflect  some  significant  features  of  the  flow,  such  as  shocks 
and  shear  layers.  It  was  very  difficult  and  time  consuming  to  construct  an  appropriate  grid.  In  partic¬ 
ular  the  feeding  sheet  was  not  properly  resolved  by  any  of  the  grids  used.  Also  very  small  time  steps, 
significant  values  of  artificial  viscosity,  and  large  numbers  of  iterations  were  required.  A  solution- 
adaptive  grid  procedure  is  developed  to  address  these  issues. 

The  coordinates  chosen  for  the  mathematical  formulation  of  fluid  flow  problems  can  have  a  signifi¬ 
cant  effect  on  the  accuracy  and  computational  requirements  of  the  analysis.  In  order  to  avoid  prob¬ 
lem-specific  flow  solvers,  it  is  necessary  that  the  coordinates  be  boundary  conforming.  Also,  it  is 
desirable  that  the  coordinate  point  distribution  reflect  some  significant  features  of  the  flow.  Typical¬ 
ly,  grids  arc  constructed  based  on  the  experience  and  intuition  of  the  user.  This  process  requires  an 
experienced  user  and  results  in  a  time-consuming,  iterative  procedure.  Several  papers  such  as  those 
by  Thompson  [1993],  Anderson  [1987],  Soni  and  Yang  [1992]  and  Soni,  Weatherill  and  Thompson 
[1993]  present  applications  where  significant  improvements  in  accuracy  have  been  obtained 
through  the  use  of  solution  adaptive  grid  procedures.  Several  approaches  have  employed  for  both 
structured  and  unstructured  grid  adaption.  The  most  widely  used  approaches  involve  grid  point  re¬ 
distribution,  local  grid  point  enrichment/derefinement  or  local  modification  of  the  actual  flow  solv¬ 
er.  However,  the  success  of  any  one  of  these  methods  ultimately  depends  on  the  feature  detection 
algorithm  used  to  flag  regions  of  the  solution  domain  which  require  a  fine  mesh  for  accurate  repre¬ 
sentation.  Typically,  weight  functions  comprised  of  combinations  of  flow  variable  derivatives  are 
constructed  to  mimic  the  local  truncation  error.  The  selection  of  which  flow  variables  and  in  what 
combinations  to  include  in  such  weight  functions,  may  require  substantial  user  input.  These  weight 
functions  can  then  be  used  to  construct  blending  functions  for  algebraic  redistribution,  interpolation 
functions  for  unstructured  grid  generation,  forcing  functions  to  attract/repel  points  in  an  elliptic  sys¬ 
tem,  or  to  trigger  local  refinement.  Dwyer  [1984]  and  Soni  and  Yang  [1992]  have  developed  scaling 
procedures  to  lessen  the  required  user  input  for  construction  of  these  weight  functions.  It  is  believed 
that  an  adaptive  grid  system  which  requires  no  user  input  while  ensuring  an  efficient  grid  point  dis¬ 
tribution  would  dramatically  increase  the  routine  use  of  CFD  in  the  design  and  analysis  environment. 
The  development  of  such  a  system  is  the  ultimate  objective  of  this  work. 

Use  of  the  adaptive  process  with  the  NEARC  flow  solver  did  improve  the  convergence  level,  agree¬ 
ment  with  flow  visualization,  and  the  agreement  with  the  measured  force  and  moment  coefficients. 
The  promising  initial  success  indicates  that  the  enhancements  proposed  significant  further  improve¬ 
ments  can  be  obtained.  Many  problems  of  engineering  interest  involve  multi-block  grids.  Hence, 
it  is  advantageous  that  the  adaptive  grid  system  be  developed  to  recognize  flow  structures  of  differ¬ 
ent  types  as  well  as  differing  intensity,  and  adequately  address  scaling  and  normalization  across 
blocks.  The  accurate  prediction  of  unsteady  flow  phenomena  is  becoming  increasingly  important. 
It  would  be  desirable  to  utilize  the  benefits  of  solution  adaptive  gridding  for  simulation  of  unsteady 
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flow  problems.  These  benefits  are  largest  when  the  physical  locations  of  the  flowfield  which  require 
a  fine  grid  for  adequate  resolution  vary  greatly  wi&  time. 

2.  DESCRIPTION  OF  PROBLEM 

In  order  to  study  the  accuracy  and  efficiency  improvements  due  to  the  grid  adaptation,  it  is  necessary 
to  quantify  grid  size  and  distribution  requirements  as  weU  as  computational  times  of  non-adapted 
solutions.  Hence,  the  first  phase  of  this  effort  involves  benchmarking  the  predictive  capability 
through  the  use  of  solution  adaptive  gridding.  This  adaptive  grid  capability  should  automatically 
resolve  complex  flows  with  shock  waves,  expansion  waves,  shear  layers  and  complex  vortex-vortex 
and  vortex-surface  interactions.  An  adaptive  grid  approach  seems  well  suited  for  such  problems 
in  which  the  spatial  distribution  of  length  scales  is  not  known  a  priori. 

There  is  interest  in  prediction  capabilities  for  vortical  flows  about  slender  bodies  of  revolution.  Ma¬ 
neuvering  missiles  often  operate  under  flight  conditions  where  forces  induced  by  vortices  often 
dominate.  Also,  delta  wings  proposed  for  High  Speed  Civil  Transport  (HSCT)  aircraft  often  utilize 
vortex  lift  during  landing  and  takeoff.  Hence,  the  design  of  such  vehicles  would  greatly  benefit  from 
a  detailed  and  accurate  knowledge  of  such  flow  fields.  Due  to  the  current  interest  in  this  type  of  flow, 
a  generic  missile  body  has  been  chosen  as  a  test  problem.  As  angle  of  attack  and  supersonic  free- 
stream  conditions,  this  configuration  exhibits  large  scale  separated  vortical  flow,  vortex-vortex  and 
vortex-surface  interactions,  separated  shear  layers  and  multiple  shocks  of  different  intensity.  Also, 
there  is  a  substantial  amount  of  experimental  data  available  over  a  wide  range  of  flow  parameters. 
In  fact,  this  problem  is  the  focus  of  the  KTA2-1 2  program.  Pagan  and  Molton  [1991],  Hsieh,  Ward- 
law  and  Birch  [1991],  as  weU  as  the  Army  Research  Laboratory  (KTA2-12)  have  compiled  experi¬ 
mental  data.  Data  for  the  Mach  number  ranges  from  0.7  to  3.5  and  angles  of  attack  from  0  to  20 
degrees  are  available.  Surface  pressure  coefficients  as  weU  as  pitot  pressure  data  for  the  off  surface 
region  are  available.  The  availability  of  detailed  experimental  data  for  verification  of  the  results 
coupled  with  the  wide  range  of  flow  features  provides  a  critical  test  for  the  capabilities  of  the  flow 
solver  used  and  any  associated  adaptive  grid  technique.  Figure  1  iUustrates  the  geometric  configura¬ 
tion  chosen  for  study  along  with  an  associated  grid  constructed  using  the  GENIE++  [Soni,  Thomp¬ 
son  et  al.,  1992]  grid  generation  system. 

3.  RESEARCH  OBJECTIVES  &  APPROACH 

An  overaU  objective  of  this  work  is  to  perform  simulations  in  support  of  the  KTA2-12  program  and 
to  address  research  issues  associated  with  grid  adaptation  arising  from  these  simulations. 

The  following  tasks  are  performed  to  meet  the  aforestated  objectives;: 

The  first  task  involves  simulations  in  order  to  better  quantify  the  predictive  capability  of  the  NPARC 
[NASA  LeRC  1993]  and  CFL3D  [NASA  LaRC  1993]  flow  solvers.  The  model  configuration  for 
the  KTA2-12  project  is  a  3-caliber  ogive  with  a  10-caliber  cylindrical  afterbody.  The  experimental 
data  is  available  for  six  test  cases  involving  different  Mach  numbers  (in  the  range  of  0.7  to  3.5)  and 
angle  of  attack  (8  and  14).  The  Navier-Stokes  simulations  are  performed  with  different  turbulence 
models.  The  six  cases  for  which  wind  tunnel  test  conditions  are  available  are  presented  in  Table  1. 

The  second  task  is  to  develop  an  adaptive  grid  strategy  and  perform  simulations  using  grid  adapta¬ 
tion.  The  results  of  grid  adaptation  should  be  compared  for  improvements  with  the  experimental 
data. 
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Case  Number 

Mach  Number 

Angle  of  Attack 

Reynolds  No,  x  10^ 

/Ft 

1 

1.45 

14 

2.0 

2 

1.8 

14 

2.0 

3 

2.5 

14 

4.0 

4 

3.5 

8 

4.0 

5 

3.5 

14 

2.0 

6 

0.7 

14 

2.0 

Table  1.  Wind  Tunnel  Test  Conditions 

4.  SIMULATION  RESULTS 

The  NPARC  system  simulates  Reynolds  averaged  Navier-Stokes  equations  using  finite  difference 
Steger-Warming  scheme  utilizing  multiblock  structured  grids.  The  scheme  is  second  order  accurate 
with  ADI  based  time  differencing  allowing  local  time  stepping.  The  program  facilitates  user  ori¬ 
ented  boundary  condition  specification.  The  CFL3D  system  also  simulates  Reynolds  averaged  Nav¬ 
ier-Stokes  equations  using  finite  volume  algorithm.  ITie  program  allows  a  choice  of  Steger-Warm¬ 
ing  or  Roe  upwinding  scheme.  The  third  order  accuracy  can  be  achieved  using  MUSCL  scheme. 
The  limiters  are  based  on  the  minmod  technique.  The  program  facilitates  local  time  stepping,  mesh 
sequencing  and  multigrid  options  to  accelerate  convergence.  The  Balwin-Lomax,  Bdwin-Barth 
and  A:— 0)  turbulence  models  are  applicable  for  simulation. 

The  solutions  were  obtained  using  both  NPARC  [NASA  LeRC  1994]  and  CFL3D  [NASA  LaRC 
1993].  Four  turbulence  models  were  used  to  study  the  effect  of  the  turbulence  model  on  the  solu¬ 
tions.  The  models  used  where  the  Baldwin-Lomax  with  the  Dgani-Schiff  correction,  Baldwin- 
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Barth,  k-2,  k—(a  models  for  NPARC  and  Baldwin-Barth,  Baldwin-Lomax,  k—o)  and  Spalart-All- 
maras  models  were  used  forCFL3D.  Significant  variations  in  the  off  surface  structure  were  apparent 
with  all  flow  solvers.  The  vortex  grows  in  strength  with  the  streamwise  coordinate,  and  is  located 
farther  from  the  surface.  The  boundary  layer  separates  to  form  the  primary  and  secondary  vortices. 
Hence,  as  the  distance  of  the  shear  layer  from  the  surface  increases,  the  grid  resolution  decreases 
and  accuracy  deteriorates  in  the  downstream  direction.  It  was  found  that  the  solutions  were  very 
sensitive  to  the  turbulence  model.  The  general  comments  on  both  these  codes  is  provided  in  Table 
2. 

The  L2  residual  is  considered  for  checking  the  order  of  the  convergence.  Most  of  the  cases  NPARC 
takes  6000-7000  iterations  to  achieve  a  residual  of  10^-10“^.  For  most  of  the  turbulence  model 
studies  the  NPARC  code  is  restarted  from  the  previously  attained  convergence  solution  for  other 
models.  This  procedure  reduces  the  overall  run  time  for  different  turbulence  model  studies. 

Typical  convergence  histories  for  different  cases  are  shown  in  Figures  2-5. 

The  results  obtained  by  applying  various  turbulence  models  for  case  3  studies  are  presented.  Four 
turbulence  models  were  used  for  NPARC  and  CFL3D  codes,  as  explained  in  the  code  performance 
section.  k-O)  implementation  in  CFL3D  code  may  be  corrupted  due  to  errors  (programming  errors). 
The  surface  presents  results  associated  with  different  turbulence  models  and  comparison  with  exper¬ 
imental  data  for  both  the  codes  are  displayed  in  Figures  6a-b  and  7a-b. 

For  CFL3D,  Spalart-Allmaras  model  gave  good  results.  Comparison  of  resulting  simulation  using 
CFL3D  and  NPARC  is  also  preformed.  These  results  are  presented  in  Figure  8a-b  for  Baldwin- 
Barth  turbulence  model,  in  Figure  8a-b  for  Baldwin-Lomax  turbulence  model  and  in  Figure  lOa-b 
for  k— 0)  model. 

The  pitot-pressure  with  vortex  structure  comparison  with  experimental  data  is  demonstrated  in  Fig¬ 
ures  11-15  for  different  axial  locations. 

The  forces  and  moments  are  calculated  for  case  3  and  is  given  in  Table  3.  The  computed  axial  forces 
are  not  matching  well  with  that  of  experiment  data.  But  simulated  the  normal  force  and  moment 
matches  well  with  the  experimental  data.  The  NPARC  code  with  Baldwin-Barth  turbulence  model 
predicted  the  normal  force  and  moments  better  than  the  other  combinations. 
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NPARC 


CFL3D 


Central  difference  with  artificial  dissipation 

Sensitive  to  grid  resolution 
esp  near  shock 

Vortex  feeding  sheet  difficult  to  resolve  properly 
Adaptive  grid  affects  solution 

Very  small  CFL  limit  (0.5) 

Apparently  due  to  shock  resolution 

Slow  convergence  (5000  iter) 
esp  in  low  speed  viscous  regions 

Cray  C90, 1,556,280  grid  points  Case3  BL 
47,332,701  words 
5894  CPU  seconds 
1.2  X  lU^  sec/grid  point 
30  words/grid  point 
275  MFLOPS 


Flux-Difference  split,  upwind 
Quick  and  easy 
Sharp  resolution  of  shocks 
Sensitive  to  Flux  limiter  used 
Not  as  sensitive  to  grid  as  NPARC 

Cray  C90, 1,556,280  grid  points  Case3  BL 
91386581  words 
5808  CPU  seconds 
9.33  X 10  sec/grid  point 
59  words/grid  point 

288  MFLOPS 


Table  2.  General  Comments  on  Flow  Solvers 


8 


Convergence  History 

Case  2 


Convergence  History  gor  Case  2  Using  NPARC  ( Laminar ) 


Convergence  History 

Case  3 
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Figure  3  Convergence  History  for  Case  3  Using  NPARC  ( BL  Model ) 


Convergence  History 

Case  6 


Convergence  History  for  Case  6  Using  NPARC  ( Euler ) 


Convergence  History 

Cases  using  CFL3D 


Convergnce  History  for  Case  3  Using  CFL3D 


Forces 

Experimental 

Data 

NPARC 

(BL) 

NPARC 

(BB) 

CFL3D 

BL 

CFL3D 

BB 

Axial 

0.1957 

0.3307 

0.3297 

0.3314 

0.3311 

Normal 

1.9100 

1.8855 

1.9098 

1.8729 

1.9148 

Moment 

10.2417 

10.0314 

10.2666 

9.9127 

10.3315 

Table  3.  Forces/Moments  ( Case  3 ) 
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‘e  8(a)  Baldwin-Barth  TM  for  Both  Codes 


XID  =  2.4 


FGigure  9(a)  Baldwin-Lomax  TM  Showing  Comparison  Between  Both  Codes 


CFL3D 


Baldwin-Lomax  TM  showing  Comparison  Between  Both  Codes 


Comparison  Between  Both  Codes 


CFL3D 


Comparison  Between  Both  Codes 


Experimental 


Computational 


Figure  11  Comparison  of  Experimental  and  Computational 
Vortex  Core  for  Case  3  at  x/d  =  5.5 


Experimental 


Computational 


Figure  12  Comparison  of  Experimental  and  Computational 
Vortex  Core  for  Case  3  at  x/d  =  11.5 


Experimental  Computational 


Experimental  Computational 


Figure  14  Comparison  of  Experimental  and  Computational 
Vortex  Core  for  Case  2  at  x/d  =  8.5 
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Experimental  Computational 
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5.  GRID  ADAPTATION 


An  adaptive  grid  system  is  developed  to  automatically  redistribute  grid  lines  based  on  an  evolving 
flow  solution.  Solution  adaptive  schemes  have  demonstrated  the  ability  to  sharply  capture  shocks. 
Hence,  in  this  work  primary  emphasis  has  been  placed  upon  the  ability  of  the  procedure  to  capture 
in  detail  the  feeding  sheet  and  vortex  definition.  Using  traditional  grid  generation  methods  it  is  diffi¬ 
cult  and  very  time  consuming  to  align  and  cluster  the  grid  with  the  feeding  sheet,  even  if  the  location 
could  be  determined  a  priori.  Appropriate  clustering  for  crisp  definition  of  the  vortex  core  is  also 
difficult.  It  is  believed  that  if  the  available  grid  points  were  distributed  in  a  near  optimum  manner 
then  all  seven  grids  would  provide  adequate  spatial  resolution.  Use  of  the  adaptive  procedure  has 
shown  substantial  improvement  in  flow  solver  convergence  behavior  and  agreement  with  flow  visu¬ 
alization  pictures. 


During  this  study,  an  adaptive  grid  system  capable  of  automatically  resolving  complex  flows  with 
shock  waves,  expansion  waves,  shear  layers  and  complex  vortex-vortex  and  vortex-surface  interac¬ 
tions  has  been  developed.  The  weight  functions  developed  utilize  scaled  derivatives  and  normaliz¬ 
ing  procedures  to  minimize  or  eliminate  the  need  for  user  input.  The  grid  redistribution  scheme  is 
based  on  the  elliptic  generation  system  with  control  functions  governed  by  the  assignment  function. 

The  elliptic  generation  system: 

=  0  (1) 

i=l  j=l  1=1 

where  r  :  Position  vector, 

gy  :  Contravariant  metric  tensor 
:  Curvilinear  coordinate,  and 
Pjfc  :  Control  function. 


is  widely  used  for  grid  generation  [Thompson  and  Warsi  1985] .  Control  of  the  distribution  and  char¬ 
acteristics  of  a  grid  system  can  be  achieved  by  varying  the  values  of  the  control  functions  Pk  in  Equa¬ 
tion  1.  The  application  of  the  one  dimensional  form  of  Equation  1  combined  with  equidistribution 
of  the  weight  function  results  in  the  definition  of  a  set  of  control  functions  for  three-dimensions. 


(/=  1,2,3) 


(2) 


These  control  functions  were  generalized  by  Eiseman  [1983]  as: 


O' =  1.2,3) 


Wi 


(3) 


In  order  to  conserve  some  of  the  geometrical  characteristics  of  the  original  grid  (he  definition  of  the 
control  functions  is  extended  as: 


Pi  {Pinilialgeometry)  +  Ci(Pj  a  =  1,2,3) 


(4) 


where  P  initial  geometry  ■  Control  function  based  on  initial 

geometry 

Pwt  :  Control  function  based  on 
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Ci 


current  solution 
:  Constant  weight  factor. 


These  control  functions  are  evaluated  based  on  the  current  grid  at  a  given  time-step.  This  can  be 
formulated  as: 

p(n)  ^  p(n-l)  +  (i  =  1,2,3)  (5) 

where 

/,(!)  =  p(o)  +  (i  =  1,2,3)  (6) 


A  flow  solution  is  first  obtained  with  an  initial  grid.  Then  the  control  function  Pj  is  evaluated  in 
accordance  with  Equations  2  and  5,  which  is  a  combination  of  the  geometry  of  the  current  grid  and 
the  weight  functions  associated  with  the  current  flow  solution. 

The  software  in  its  current  form  inputs  two  PLOT3D  format  files  [Buning  1985],  one  for  the  grid 
and  one  for  the  flow  solution.  Output  consists  of  an  NPARC  restart  file,  as  well  as  two  PLOT3D 
files  of  the  adapted  grid  and  the  flow  solution  interpolated  onto  the  new  grid.  A  multi-block  treat¬ 
ment  capability  is  included,  but  has  not  yet  been  validated.  The  adaptive  grid  is  constructed  in  three 
steps.  The  first  step  is  to  generate  the  weight  functions,  which  due  to  their  critical  importance  will 
be  discussed  in  detail  in  a  separate  section.  The  second  step  is  to  generate  the  actual  adapted  grid 
by  equidistribution  of  the  aforementioned  weight  function.  In  the  current  work  this  is  accomplished 
by  the  numerical  solution  of  Equation  1 .  A  coupled  three-dimensional  strongly  implicit  procedure 
(CSEP)  as  described  by  Ghia  et  al.  [1981]  has  been  implemented  for  the  solution  of  the  discretized 
equations.  Upwind  differencing,  with  biasing  based  on  the  sign  of  the  forcing  functions,  as  well  as 
central  differencing  has  been  implemented  and  studied  for  the  first  derivative  terms.  The  first  order 
upwind  differencing  increases  the  stability  of  the  procedure,  but  at  the  expense  of  smearing  the  grid 
clustering.  Hence,  a  hybrid  upwind/central  differencing  scheme  has  been  implemented  to  lessen  this 
smearing  while  maintaining  the  smoothness  and  stability  of  the  upwind  procedure.  Central  differ¬ 
encing  has  been  employed  for  all  second  and  mixed  derivative  terms.  AH  non-linear  terms  were 
treated  by  quasi-Unearization.  Boundary  point  movement  was  allowed  through  Neumann  boundary 
conditions  for  plane  surfaces,  and  by  a  NURBS  fit  and  redistribution  for  curved  surfaces.  The 
NURBS  redistribution  is  adapted  from  the  GENIE++  code.  For  the  curved  surfaces  the  redistribu¬ 
tion  mesh  was  based  on  the  nearest  interior  coordinate  surface.  Since  Equation  1  is  solved  iterative¬ 
ly,  the  forcing  functions  must  be  evaluated  for  each  new  successive  grid  point  location.  The  forcing 
functions  are  obtained  by  interpolation  from  the  original  grid.  The  interpolation  procedure 
employed  was  that  used  for  the  non-matching  block  to  block  interface  capability  of  the  NPARC 
code.  Upon  convergence  of  this  process  the  flow  solution  is  interpolated  onto  the  adapted  grid  using 
the  same  interpolation  subroutine.  Calculations  can  then  be  continued  from  the  new  restart  file .  This 
procedure  can  then  be  repeated  until  an  acceptable  solution  is  obtained.  Experience  indicates  that 
coupling  this  procedure  with  a  code  capable  of  treating  time  accurate  grid  movement  would  ease 
this  process  and  lessen  the  CPU  requirements. 

Application  of  the  equidistribution  law  results  in  grid  spacing  inversely  proportional  to  the  weight 
function,  and  hence,  the  weight  function  determines  the  grid  point  distribution.  Ideally,  the  weight 
would  be  the  local  truncation  error  ensuring  a  uniform  distribution  of  error.  For  a  given  discretiza¬ 
tion  the  truncation  error  term  contains  derivatives  of  order  greater  than  that  of  the  solution  order  of 
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accuracy.  Evaluation  of  higher-order  derivatives  from  discrete  data  is  progressively  less  accurate 
and  subject  to  noise.  However  lower-order  derivatives  must  be  non-zero  in  regions  of  wide  varia¬ 
tion  of  higher-order  derivatives,  and  are  proportional  to  the  rate  of  variation.  Therefore,  it  is  pos¬ 
sible  to  employ  lower-order  derivatives  as  a  proxy  for  the  truncation  error.  Determination  of  this 
function  is  one  of  the  most  challenging  areas  of  adaptive  grid  generation. 


Analysis  of  the  weight  functions  explored  to  date  indicates  that  density  or  velocity  derivatives  inde¬ 
pendently  are  not  sufficient  to  represent  the  different  types  and  strengths  of  flow  features  for  viscous 
flows.  Density  or  pressure  for  that  manner  varies  insufficiently  in  the  boundary  layer  to  be  used  to 
construct  weight  functions  for  representation  of  these  features.  While  velocity  derivatives  by  them¬ 
selves,  for  viscous  flows,  are  dominated  by  the  boundary  layer,  and  additional  variables  must  be  in¬ 
cluded  to  represent  other  flow  features.  The  present  weight  function  consists  of  relative  derivatives 
of  density,  and  the  three  conservative  velocities. 


= 


Qijt 


(7) 


where  q  =  {p,pu,£>v,pw} 


The  relative  derivatives  are  necessary  to  detect  features  of  varying  intensity,  so  that  weaker,  but  im¬ 
portant  structures  such  as  vortices  are  accurately  reflected  in  the  weight  function.  One-sided  differ¬ 
ences  are  used  at  boundaries,  and  no-slip  boundaries  require  special  treatment  since  the  velocity  is 
zero.  This  case  is  handled  in  the  same  manner  as  zero  velocity  regions  in  the  field.  A  small  value, 
epsilon  in  the  preceding  equation,  is  added  to  all  normalizing  quantities.  Also  it  appears  that  the 
Boolean  sum  construction  method  of  Soni  and  Yang  [1992]  would  more  evenly  balance  the  weight 
functions,  as  several  features  are  reflected  in  multiple  variables,  while  some  are  reflected  in  only  one. 


The  weight  function  developed  to  represent  the  error  associated  with  the  discretization  is  presented 
in  Figures  16, 17,  and  18  for  the  NPARC,  CFL3D  and  adapted  NPARC  calculations  respectively. 
As  can  be  seen  from  Figures  16, 17,  and  18  the  feeding  sheet,  and  both  primary  and  secondary  vor¬ 
tices,  as  well  as  the  strong  shock  at  the  nose  are  clearly  visible.  It  can  be  seen  in  Figure.  18  that  the 
solution  has  become  much  smoother  and  both  the  shock  and  feeding  sheet  are  much  sharper  and 
clearer.  Only  structures  that  have  been  at  least  partially  resolved  by  the  flow  solver  can  be  detected 
by  the  weight  function.  Hence  the  quality  of  the  weight  function  is  dependent  upon  the  quality  of 
the  solution. 

Flow  at  Mach  number  1 .45  and  14  degree  angle  of  attack  has  been  simulated  around  the  missile  using 
the  NPARC  flow  solver.  The  grid  constructed  for  the  second  cycle  of  adaption  using  the  hybrid  dif¬ 
ferencing  scheme  for  the  grid  equations  is  presented  in  Figure  19.  Figure  20  presents  a  side  by  side 
comparison  of  the  grid  after  two  and  three  adaptation  cycles.  The  alignment  of  the  grid  lines  with 
the  flow  structures  of  interest  is  clearly  visible.  It  can  be  seen  that  the  grid,  and,  hence,  the  solution 
varies  substantially  from  the  second  to  the  third  cycle.  The  adaption  procedure  and  the  flow  solver 
should  be  coupled  so  that  the  evolving  grid  can  continuously  reflect  features  that  are  detected  as  the 
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solution  develops  and  improves  due  to  adaptation.  This  should  increase  the  efficiency  of  the  proce¬ 
dure. 


Figure  16  Weight  Function  Based  Upon  Original  NPARC  Solution. 


Figure  17  Weight  Function  Based  Upon  Original  CFL3D  Solution. 

The  results  shown  demonstrate  the  capability  of  the  procedure  to  detect  and  resolve  shocks  of  differ¬ 
ing  strengths,  primary  and  secondary  vortices,  and  shear  layers  adequately.  The  use  of  adaptive  re- 
meshing  allows  the  use  of  larger  time  steps  resulting  in  increased  convergence.  However,  the  single 
greatest  benefit  resulting  from  the  adaptive  procedure  is  a  dramatic  lowering  of  the  artificial  viscos¬ 
ity  parameters  required  for  stability.  This  results  in  off  surface  structures  which  are  much  sharper 
and  more  closely  resemble  the  flow  visualization.  No  special  treatment  is  needed  for  boundary  layers 
and  no  user  input  is  required.  In  fact,  the  normal  spacing  for  the  first  point  off  the  wall  on  the  wind¬ 
ward  side  was  reduced  to  on  the  order  of  10~^  from  on  the  order  of  1(H  by  the  adaptation  procedure. 
It  should  be  noted  that  this  is  due  primarily  to  the  decrease  in  boundary  layer  thickness  as  the  artificial 
dissipation  values  were  lowered. 

As  discussed  previously,  it  is  imperative  that  the  adaptation  process  be  coupled  with  the  flow  solver 
as  experience  indicates  that  complex  flowfields  may  require  dozens  of  adaptive  cycles.  This  is  par- 


32 


Figure  20  Comparison  of  adapted  grid  after  two  and  three  cycles. 


ticularly  critical  for  flow  fields  that  are  initially  poorly  resolved.  Also,  the  interpolation  procedure 
for  updating  the  forcing  functions,  which  is  very  time  consuming  and  often  a  source  of  trouble, 
would  be  eliminated. 

Currently,  simple  multi-block  problems  are  being  simulated  to  evaluate  this  capability.  Attention 
is  being  placed  upon  across  block  scaling  for  the  weight  functions.  A  global  maximum  across  all 
blocks  seems  to  work  well  for  the  normalization  of  weight  function  components.  However,  this  is¬ 
sue  is  being  monitored  as  more  complex  problems  are  attempted.  A  further  problem  arises  due  to 
equidistribution  via  forcing  functions.  This  does  not  appear  to  be  a  problem  where  a  block  face  is 
connected  to  one  and  only  one  face  of  another  block.  Problems  have  been  encountered  for  a  multi¬ 
block  launch  vehicle  computation  where  a  block  face  consisted  of  a  block  to  block  connection  and 
a  solid  surface  for  the  base  region.  A  smooth,  continuous  set  of  forcing  functions  across  this  block 
boundary  does  not  yield  a  smooth  grid.  It  is  believed  that  this  is  due  to  the  elliptic  character  of  the 
grid  equations.  The  simple  fix  is  to  introduce  artificial  block  boundaries  to  force  one  to  one  corre¬ 
spondence. 

Use  of  the  adaptive  process  with  the  NPARC  flow  solver  did  improve  the  convergence  level,  agree¬ 
ment  with  flow  visualization,  and  the  agreement  of  the  force  and  moment  coefficients  with  the  ex¬ 
perimental  data.  The  promising  initial  success  indicates  that  with  the  enhancements  proposed  signif¬ 
icant  further  improvements  can  be  obtained.  It  is  proposed  to  work  on  block  interface  issues  as  well 
as  coupling  with  flow  solvers,  unsteady  flows,  and  more  complex  multi-block  problems. 

The  simulation  results  for  grid  adaptation  are  presented  in  Figures  21-31. 

From  the  grid  adaption  study  it  was  concluded  that  the  results  are  different  than  the  unadapted  one. 
However,  not  clear  improvement  in  matching  experimental  pressure  data.  The  vortex  sheet  and  the 
shock  structure  is  better  pronounced  in  adaptive  solution  is  in  good  agreement  with  the  experimental 
data.  Further  work  is  needed  in  grid  adaptation  and  grid  optimization.  However,  there  is  a  remark¬ 
able  improvement  in  the  simulated  values  of  normal  force  and  moment  coefficients  using  adapted 
grid  as  compared  to  the  unadapted  grid.  The  calculated  values  of  these  coefficients  are  compared 
with  the  experimental  data  in  Table  4. 


Forces 

Experimental 

Data 

Unadapted 
(NPARC,  BL) 

Adapted 
(NPARC,  BL) 

Axial 

0.1957 

0.3307 

0.3309 

Normal 

1.9100 

1.8855 

1.9052 

Moment 

10.2417 

10.0314 

10.2666 

Table  4.  Forces/Moments  Comparison 
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Figure  21  Comparison  of  Unadapted  and  Adapted  Grid 


Unadapted  H  0.2  r  -  Unadapted 


Comparison  of  Unadapted  and  Adapted  Solution  for  Case  3 


Unadapted  A  0,1  k*.  -  Unadapted 
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Figiire  23(a)  Comparison  of  Unadapted  and  Adapted  Solutions  for  Case  2 


Unadapted  -  -  Unadapted 

Adq?ted  0.1  -  . Adapted 


■e  23(b)  Comparison  of  Unadpted  and  Adapted  Solutions  for  Case  2 


ST  -  a/X  OHVdM  9  ^^0  P'Z  -  Q/X  DUVdN  9 


40 


Figure  24(a)  Comparison  of  Unadapted  and  Adapted  Solution  for  Case  6  ( Euler ) 


Unadapled  -  Unadapted 

Adapted  -  0,1  -  . Adapted 
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Figure  24(b)  Comparison  of  Unadapted  and  Adapted  Solutions  for  Case  6  ( Euler ) 


Figure  25  Comparison  of  Experimental  and  Computational 
(Adapted)  Vortex  Core  for  Case  3  at  x/d  =  5.5 
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Unadapted  Adapted 


Figure  26  Comparison  of  Unadapted  and  Adapted 
Vortex  Core  for  Case  3  at  x/d  =  5.5 


Figure  28  Comparison  of  Unadapted  and  Adapted  Vortex 
Core  for  Case  3  at  x/d  =  11.5 
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Figure  29  Comparison  of  Unadapted  and  Adapted 
Vortex  Core  for  Case  3  at  x/d  =  5.5 


Figure  30  Comparison  of  Unadapted  and  Adapted 
Vortex  Core  for  Case  3  at  x/d  =  5.5 


Figure  31  Comparison  of  Unadapted  and  Adapted 
Vortex  Core  for  Case  3  at  x/d  =  5.5 


6.  CONCLUSION 


The  CFD  simulations  performed  utilizing  the  flow  solvers:  NPARC  and  CFL3D  were  submitted 
to  the  KTA2-12  team.  Nine  different  CFD  codes  were  exercised  by  the  KTA  team  in  all  six  cases. 
However,  the  concentration  was  placed  on  case  3  (Mach=2.5,  AOA=14)  for  the  reporting  and  com¬ 
parison  purpose.  Qualitatively,  consistent  results  were  obtained  by  all  these  codes.  Also,  no  best 
turbulence  model  could  be  identified  as  of  result  of  these  simulations  in  view  of  comparison  with 
experimental  data.  A  technical  paper  summarizing  this  program  has  been  written  and  was  presented 
at  the  AIAA  conference.  The  paper  entitled,  “The  Application  of  CFD  to  the  Prediction  of  Missile 
Body  Vortices”  is  included  in  the  Appendix  A.1  of  this  report. 

An  adaptive  grid  strategy  has  been  developed  and  grid  adaptation  was  performed  to  evaluate  the  in¬ 
fluence  on  simulation.  Little  variation  was  observed  in  the  surface  pressure  coefficients  due  to 
adaptation.  However,  the  vortex  core  resolution  was  crisp  and  the  forces  and  moments  indicated 
better  agreement  because  of  adaptation.  Further  work  in  grid  adaptation  is  needed  to  evaluate  effect 
of  the  grid  adaptation  on  simulations.  The  development  of  weight  function  and  detailed  study  of 
grid  adaptation  has  been  reported  in  a  technical  paper  entitled,  “A  Structured  Grid  Based  Solution- 
Adaptive  Technique  for  Complex  Separated  Flow”.  This  paper  has  been  accepted  for  the  publication 
in  the  Journal  of  Applied  Mathematics  and  Computation  and  is  included  in  the  Appendix  A.2  of  this 
report. 


49 


REFERENCES 


1.  Anderson,  D.A.,  (1987),  Equidistribution  schemes,  Poisson  generators,  and  adaptive  grids, 
Appl.  Math.  Comput  24,  211-227. 

2.  Anderson,  D. A.,  “Adaptive  Grid  Methods  for  Partial  Differential  Equations,”  Advances  in 
Grid  Generation,  ASME  Fluids  Engineering  Conference,  Houston,  1983. 

3.  Benson,  Rusty  A.,  and  McRae,  D.  Scott,  “Time-Accurate  Simulations  of  Unsteady  Flows 
With  A  Dynamic  Solution-Adaptive  Grid  Algorithm,”  Numerical  Grid  Generation  in  Com¬ 
putational  Fluid  Dynamics  and  Related  Fields,  (1994)  (Pineridge  Press  Limited,  Swansea, 
U.K.)  Ed.  N.P.  Weatherill,  P.R.  Eiseman,  J.  Hauser,  J.F.  Thompson. 

4.  BockeUe,  M.J.  (1988)  Adaptive  grid  movement  schemes  and  the  numerical  simulation  of 
shock  vortex  interaction,  Ph.D.  Thesis,  Columbia  University,  New  York. 

5.  Brackbill,  J.S.,  “Application  and  Generalization  of  Variational  Methods  for  Generating 
Adaptive  MEshes,”  Numerical  Grid  Generation,  Edited  by  J.G.  Thompson,  Elsevier  Sci¬ 
ence  Publishing  Company,  1985. 

6.  Brackbill,  J.U.  and  J.  Saltzman  (1982)  Adaptive  zoning  for  singular  problems  in  two  dimen¬ 
sions,  J.  Comput.  Phys.  46,  342-368. 

7.  Buning,  P.G.  and  Steger,  J.L.,  ’’Graphics  and  Flow  Visualization  in  Computational  Fluid  Dy¬ 
namics,”  AIAA  Paper  85-1507-CP,  Proceedings  of  the  AI AA  7^  Computational  Fluid  Dy¬ 
namics  Conference.  1985. 

8.  Catherall,  D.,  “An  Assessment  of  the  Advantages  of  Adaptivity  on  Structured  Grids,”  Nu¬ 
merical  Grid  Generation  in  Computational  Fluid  Dynamics  and  Related  Fields,  (1994) 
(Pineridge  Press  Limited,  Swansea,  U.K.)  Ed.  N.P.  Weatherill,  P.R.  Eiseman,  J.  Hauser,  J.F. 
Thompson. 

9.  Dwyer,  H.A.,  ’’Grid  Adaption  for  Problems  in  Fluid  Dynamics,”  AIAA  Journal,  vol.  22.  No. 
12,  pp.  1705-1712,  December  1984. 

10.  Eiseman,  P.R.,  ’’Alternating  Direction  Adaptive  Grid  Generation,”  AIAA  Paper  83-1937, 
1983. 

11.  Ghia,  K.N.,  Ghia,  U.,  Shin,  C.T.  and  Reddy,  D.R.,  ’’Multigiid  Simulation  of  Asymptotic 
Curved-Duct  Flows  Using  a  Semi-ImpUcit  Numerical  Technique,”  Computers  in  Flow  Pre¬ 
diction  and  Fluid  Dynamics  Experiments.  ASME  Publication,  New  York  1981. 

12.  Hagmeijer,  R.,  (1993)  “Anisotropic  Grid  Adaption  Based  on  Diffusion  Equations,”  M/men- 
cal  Grid  Generation  in  Computational  Fluid  Dynamics  and  Related  Fields,  (1994)  (Piner¬ 
idge  Press  Limited,  Swansea,  U.K.)  Ed.  N.P.  Weatherill,  P.R.  Eiseman,  J.  Hauser,  J.F. 
Thompson. 

13.  Hsieh,  T.,  Wardlaw,  A.B.  and  Birch,  T.J.,  ’’Vortical  Flows  about  a  Long  Ogive-Cylinder  at 
M=3.5  and  a=18®”,  AIAA-91-1808,  AIAA  22”^  Fluid  Dynamics.  Plasma  Dynamics  and 
Lasers  Conference.  Honolulu,  Hawaii,  June  24-26, 1991. 


50 


14.  Jones,  W.P.,  Menzies,  K.R.,  “An  Adaptive  Grid  Redistribution  Scheme  and  its  Application 
to  Incompressible  Recirculating  Fluid  Flow  Calculations,”  Numerical  Grid  Generation  in 
Computational  Fluid  Dynamics  and  Related  Fields,  (1994)  (Pineridge  Press  Limited,  Swan¬ 
sea,  U,K.)  Ed.  N.P.  Weatherill,  PR.  Eiseman,  J.  Hauser,  J.F.  Thompson. 

15.  Kim,  H.J.  and  J.F.  Thompson  (1990),  Three-dimensional  adaptive  grid  generation  on  a  com¬ 
posite-block  grid,  AIAA  J.  28  470-477. 

16.  NASA  LeRC  and  USAF  AEDC,  ”NPARC  1.0  User  Notes,”  June  1993. 

17.  NASA  LaRC,  ’’User  Document  for  CFL3D/CFL3DE  (Version  1.0)”,  1993. 

18.  Niederdrenk,  Peter,  “Grid  Adaptation  to  Multiple  Auto-scaled  Solution  Features,”  Numeri¬ 
cal  Grid  Generation  in  Computational  Fluid  Dynamics  and  Related  Fields,  (1994)  (Piner¬ 
idge  Press  Limited,  Swansea,  U.K.)  Ed.  N.P.  Weatherill,  PR.  Eiseman,  J,  Hauser,  J.F. 
Thompson. 

19.  Pagan,  D.  and  Molton,  R,  ’’Basic  Experiment  on  a  Supersonic  Vortex  Flow  Around  a  Missile 
Body,”  AIAA-91-0287,  29^  Aerospace  Sciences  Meeting.  Reno,  NV,  Jan.  7-10, 1991. 

20.  Shih,  M,H.,  “Towards  A  Comprehensive  Computational  Simulation  System  for  Ttirboma- 
chinery”.  Dissertation,  Mississippi  State  University,  May,  1994. 

21 .  Soni,  B.K. ,  “Algebraic  Methods  and  CAGD  techniques  in  Structured  Grid  Generation,”  Pro¬ 
ceedings  of  the  Fourth  International  Conference  of  Numerical  Grid  in  CFD,  Swansea, 
Wales,  April  1994. 

22.  Soni,  B.K.,  Weatherill,  N.P.  and  Thompson,  J.F.,  ’’Grid  Adaptive  Strategies  in  CFD,”  In¬ 
ternational  Conference  on  Hydo  Science  and  Engineering.  Washington,  D.C.,  June  7-11, 
1993. 

23.  Soni,  B.K.  and  Yang,  J.C.,  ’’General  Purpose  Adaptive  Grid  Generation  System,” 
AIAA-92-0664, 30^  Aerospace  Sciences  Meeting.  Reno,  NV,  Jan.  6-9, 1992. 

24.  Soni,  B.K.,  Thompson,  J.F.,  Stokes,  M.L.  and  Shih,  M.H.,  ’’GENIE++,  EAGLEView  and 
TIGER:  General  and  Special  Purpose  Interactive  Grid  Systems,”  AIAA-92-0071,  30^ 
Aerospace  Sciences  Meeting.  Reno,  NV,  Jan.  6-9, 1992. 

25.  Sweby,  P.K,  and  Yee,  H.C.,  “On  The  Dynamics  of  Some  Grid  Adaption  Schemes,”  Numeri¬ 
cal  Grid  Generation  in  Computational  Fluid  Dynamics  and  Related  Fields,  (1994)  (Piner¬ 
idge  Press  Limited,  Swansea,  U.K.)  Ed.  N.P.  Weatherill,  P.R.  Eiseman,  J.  Hauser,  J.F. 
Thompson. 

26.  Takahashi,  Satoshi,  Eiseman,  Peter  R.,  “Adaptive  Grid  Movement  with  Respect  to  Bound¬ 
ary  Curvature,”  Numerical  Grid  Generation  in  Computational  Fluid  Dynamics  and  Related 
Fields,  (1994)  (Pineridge  Press  Limited,  Swansea,  U.K.)  Ed.  N.P.  Weatherill,  P.R.  Eiseman, 
J.  Hauser,  J.F.  Thompson. 

27.  Thompson,  J.F.,  and  WeatheriU,  N.P.,  “Aspects  of  Numerical  Grid  Generation:  Current  Sci¬ 
ence  and  Art,”  AlAA-93_3539-CP,AppliedAerodynamics  Conference,  Monterey,  CA,  Au¬ 
gust  1993. 


51 


28.  Thompson,  J.R,  Z.U.A.  Warsi  and  C.W.  Mastin  (1985)  Numerous  Grid  Generation: 
Foundations  and  Applications  (North-Holland,  Amsterdam). 

29.  Thompson,  J.R,  (1985),  A  survey  of  dynamically-adaptive  grids  in  the  numerical  solution 
of  partial  differential  equations,  AppZ.  Num.  Math.  1,  3-27. 

30.  Thompson,  J.R,  RC.  Thames  and  C.W.  Mastin  (1974)  Automatic  numerical  generation  of 
body-fitted  curvilinear  coordinate  system  for  fields  containing  any  number  of  arbitrary  two- 
dimensional  bodies.  J.  Comp.  Phys.  15,  pp.  299-319. 

31.  Thornburg,  H.J.  and  Soni,  B.K.,  “Weight  Functions  inGrid  Adaption,”  Proceedings  of  the 
4th  International  Conference  in  Numerical  Grid  Generation  in  Computational  Fluid  Dynam¬ 
ics  and  Related  Fields  held  at  Swansea,  Wales  6-8th  April  1994. 

32.  Thornburg,  H.J.  (1991)  An  adaptive-grid  technique  for  simulation  of  complex  unsteady 
flows,  Ph.D.  Thesis,  Dept,  of  Mechanical,  Industrial  and  Nuclear  Engineering,  University 
of  Cincinnati,  Cincinnati,  Ohio. 

33.  Trompert,  R.A.  and  J.G.  Verwer  (1990),  A  static-regridding  method  for  two-dimensional 
parabolic  partial  differential  equations.  Report  NM-R8923,  Centre  for  Mathematics  and 
Computer  Science,  Amsterdam,  The  Netherlands. 

34.  Tu,  Y.  and  Thompson,  J.R,  “Three-Dimensional  Adaptive  Grid  Generation  on  Composite 
Configurations,”  AIAA /oMma/,  Vol.  29,  p.  2025,  1991. 

35.  Vatsa,  V.N.,  Sanetrik,  M.D.,  and  Parlette,  E.B.,  “Development  of  a  Flexible  and  Efficient 
Multigrid  Based  Flow  Solver,”  AIAA-93-0677,  Reno,  NV  Jan.  1993. 

36.  Yang,  J.C,,  and  Soni,  B.K.,  “Structured  Adaptive  Grid  Generation,”  Proceedings  of  the  Mis¬ 
sissippi  State  University  Annual  Conference  on  Differential  Equations  and  Computational 
Simulation,  MSU,  March  19-20, 1993. 

37.  Yu,  Tzu-Yi  and  Soni,  Bharat  K.,  “Geometry  Transformer  and  NURBS  in  grid  generation”. 
The  Fourth  International  Conference  in  Numerical  Grid  Generation  and  related  fields, 
Swansea,  England,  April  6-8, 1994. 

38.  Yu,  T.Y.,  “IGES  Transformer  and  NURBS  In  Grid  Generation,”  Master’s  Thesis,  Mississippi 
State  University,  August. 


52 


APPENDIX  A.1 

Technical  Paper  AIAA-97-0637 


53 


AIAA  97-0637 

The  Application  of  CFD  to  the  Prediction  of 
Missile  Body  Vortices 

Walter  B.  Sturek ,  U.  S.  Army  Research  Laboratory,  Aberdeen  Proving  Ground,  MD 

Trevor  Birch,  Defense  Research  Agency,  Bedford,  UK 

Marc  Lauzon,  Defense  Research  E^tablishment-Valcartier,  Valcartier,  CA 

Clinton  Housh,  U.  S.  Naval  Air  Warfare  Center,  China  Lake,  CA 

Joseph  Manter,  U.  S.  Air  Force  Wright  Laboratory,  Wright-Patterson  AFB,  OH 

Eswar  Josyula,  U.  S.  Air  Force  Wright  Laboratory,  Wright-Patterson  AFB,  OH 

Bharat  Soni,  Mississippi  State  University,  Starkville,  MS 


35th  Aerospace  Sciences 
Meeting  &  Exhibit 

January  6-10, 1997  /  Reno,  NV 


For  permission  to  copy  or  republish,  contact  the  American  Institute  of  Aeronautics  and  Astronautics 
1801  Alexander  Bell  Drive,  Suite  500,  Reston,  VA  22091 


THE  APPLICATION  OF  CFD 
TO  THE  PREDICTION  OF  MISSILE  BODY  VORTICES 

Waller  B.  Sturek* ,  U.  S.  Army  Research  Laboratory,  Aberdeen  Proving  Ground,  MD 
Trevor  Birch,  Defense  Research  Agency,  Bedford,  UK 
Marc  Lauzon,  Defense  Research  Establishment- Valcartier,  Valcartier,  CA 
Clinton  Housh,  U.  S.  Naval  Air  Warfare  Center,  China  Lake,  CA 
Joseph  Manter** ,  U.  S.  Air  Force  Wright  Laboratory,  Wright-Patterson  AFB,  OH 
Eswar  Josyula** ,  U.  S.  Air  Force  Wright  Laboratory,  Wright-Patterson  AFB,  OH 
Bharat  Soni*,  Mississippi  State  University,  Starkville,  MS 


Abstract 

This  paper  reports  the  results  of  a  collaborative  study 
carried  out  under  the  auspices  of  The  Technical 
Cooperation  Program  (TTCP)  with  participants  from 
Canada,  UK,  and  the  US.  The  purpose  of  this  study 
was  to  apply  Navier-Stokes  computational 
techniques  to  a  complex  flow  field  with  highly 
separated  flow  for  a  missile  shape  to  evaluate  the 
predictive  technology.  Nine  Navier-Stokes 
computational  codes  were  applied  to  predict  the  flow 
about  an  ogive-cylinder  configuration  for  transonic 
and  supersonic  velocities  at  8°  and  14°  angle  of 
attack.  The  computational  results  were  compared  to 
experimental  measurements  for  surface  pressure, 
pitot  surveys  of  the  outer  flow  field,  and  strain  gage 
force  measurements.  Computational  results  were 
obtained  for  nine  turbulence  models,  laminar  flow, 
and  inviscid  flow.  For  the  conditions  of  this  study, 
no  'best'  turbulence  model  could  be  identified  in  the 
comparisons  to  experiment.  Predictions  of  surface 
pressure  and  outer  vortical  flow  in  regions  of  highly 
separated  flow  indicate  further  development  of  the 
predictive  technology  is  required;  however,  the 
prediction  of  aerodynamic  lift  and  pitching  moment 
are  within  acceptable  design  requirements. 

Introduction 

Computational  Fluid  Dynamics  (CFD)  has  emerged 
as  a  critical  technology  for  the  design  and 
assessment  of  weapon  systems  in  the  future.  Such 
influencing  factors  as:  advancing  computer 
technology;  funding  shortfalls  that  limit  the 
availability  of  experimental  testing  in  the  design 
phase  of  weapon  systems  development;  and  advances 
in  computational  modeling  technology  have  served 


to  focus  emphasis  on  the  development  and  validation 
of  computational  predictive  technology.  A  previous 
TTCP  study,  which  was  completed  in  1990, 
addressed  a  wide  selection  of  test  cases  including 
base  flows,  bodies  of  revolution,  and  finned  missiles. 
Transonic  and  supersonic  velocities  for  angles  of 
attack  from  small  to  large  were  included  in  the  test 
matrix.  The  computational  results  were  dominated 
by  Euler  solutions.  One  result  of  this  effort  was  the 
identification  of  the  inability  of  Euler  solutions  to 
provide  adequate  accuracy  for  flows  with  significant 
boundary  layer  separation,  which  results  in  vortex 
dominated  flow.  Another  result  was  the  realization 
of  the  magnitude  of  effort  and  computer  resources 
required  to  perform  Navier-Stokes  modeling. 

Since  the  completion  of  this  previous  study, 
significant  advances  in  computational  modeling 
techniques  and  increased  availability  of 
supercomputer  resources  have  taken  place. 
Accordingly,  this  study  was  chartered  to  focus  effort 
to  apply  3-D  Navier-Stokes  computational 
techniques  to  predict  the  flow  fields  about  long  1/d 
bodies  at  moderate  angles  of  attack  for  transonic  and 
supersonic  velocities.  The  Defense  Research 
Agency,  UK  provided  a  high-quality  experimental 
data  base  (l)for  a  missile  shape  that  included  surface 
pressure,  flow  field  pitot  pressure  surveys,  and  force 
measurements.  Tbe  computational  results  have  been 
compared  to  these  experimental  data  along  with 
cross  comparisons  among  the  various  computational 
techniques. 
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Computational  Techniques 

The  participants  utilized  Navier-Stokes 
computational  codes  that  are  used  in  the 
performance  of  missile  aerodynamic  analysis  duties. 
The  codes  used  cover  a  wide  spectrum  of  predictive 
technology,  including  finite-difference,  finite- 
volume,  finite-element,  explicit,  implicit,  upwind, 
and  central-difference  formulations.  These  codes 
have  been  reported  elsewhere  in  the  literature  (2-8) 
in  detail  so  no  extensive  code  descriptions  will  be 
included  in  this  report.  Table  1  gives  a  listing  of  the 
code's  used  in  this  study. 

Experimental  Data 

The  experimental  data  were  provided  by  the  DRA, 
Bedford,  UK.  The  data  were  obtained  as  part  of  a 
comprehensive  experimental  study  that  extended 
over  a  period  of  several  years,  Birch  (1).  The  model 
configuration  of  interest  for  this  study  is  a  3-caliber 
ogive  with  a  10-caliber  cylindrical  afterbody.  A 
Schlieren  photograph  showing  the  model  mounted  in 
the  wind  tunnel  for  test  conditions  of  Mach  =  3.5 
and  alpha  =  8°  is  shown  in  Figure  1.  The  test 
conditions  for  the  six  test  cases  used  in  this  study  are 
listed  in  Table  2.  The  experimental  data  obtained 
include  surface  pressures,  pitot  surveys  of  the  flow 
field,  and  force  balance  measurements.  The  surface 
pressure  measurements  include  32  axial  stations  for 
41  circumferential  locations  from  the  wind-  to  the 
lee-side  at  a  4.5°  increment.  The  pitot  pressure 
surveys  provide  a  dense  array  of  measurements  at 
axial  stations  as  indicated  in  Table  3.  Strain  gage 
force  balance  measurements  are  available  as 
indicated. 

Results 

Test  Case  Overview 

Each  test  case  presents  significant  physical  modeling 
challenges.  The  following  discussion  attempts  to 
provide  an  overview  of  the  flow  fields  along  with  a 
“typical”  comparison  between  computation  and 
experiment.  Due  to  the  large  number  of 
computational  results,  space  limitations  prevent 
showing  a  comprehensive  comparison  between  the 
computational  techniques.  The  calculations  shown 
were  performed  using  the  overflow  code.  The 
turbulent  viscous  results  were  obtained  using  the 
BEDS  turbulence  model. 


Case  1.  Mach  =  1.45,  a  =  14°,  Re/d  =  0.667xl0^ 
Figure  2. 

The  laminar  results  indicate  strong  separation  on  the 
wind-side  prior  to  x/d  =  3.5,  whereas  the  turbulent 
results  and  experiment  indicate  separation  on  the 
lee-side  for  3.5<x/d<4.5.  The  lee-side  is  strongly 
separated  for  x/d>4.5.  The  results  at  x/d=6.5,7.5  and 
9.5  show  the  experimental  pressure  greater  than  the 
prediction  on  the  wind-side  where  close  agreement 
between  computation  and  experiment  is  expected. 
This  behavior  is,  perhaps,  evidence  that  a  wind-side 
shock  is  triggered  from  a  position  upstream  of 
x/d=6.5  and  results  in  a  pressure  jump  that  the 
computed  results  do  not  show.  An  explanation  for 
the  computation  not  predicting  this  behavior  is 
unknown.  It  is  noted  that  the  laminar  predictions 
also  do  not  predict  the  increase  in  surface  pressure 
on  the  wind-side  for  x/d>5.5  even  though  the 
laminar  predictions  indicate  early  boundary-layer 
separation  that  could  trigger  a  wind-side  shock. 

Case  2.  Mach  =  1.8,  a  =  14°,  Re/d  =  0.667x10®, 
Figure  3. 

The  laminar  results  indicate  strong  separation  for 
x/d=3.5,  which  is  not  supported  by  the  experimental 
data  or  the  turbulent  calculations.  The  laminar 
viscous  results  consistently  predict  early  flow  field 
separation  compared  to  the  experiment  and 
turbulent  viscous  results.  The  experiment  indicates 
the  flow  is  strongly  separated  for  x/d  >  4.5.  A 
pressure  jump  is  noted  at  x/d  =  9.5  similar  to  that 
noted  for  Case  1. 

Case  3.  Mach  =  2.5,  a  =  14  °,  Re/d  =  1.123x10®, 
Figure  4. 

The  laminar  viscous  results  indicate  early  separation 
^  of  the  viscous  layer  compared  to  the  experiment  and 
the  turbulent  viscous  results.  Good  agreement  is 
achieved  on  the  wind-side  for  all  results  shown. 

This  test  case  was  selected  as  the  top  priority  because 
the  Mach  number  and  Reynolds  number  of  the  flow 
provided  the  expectation  for  high-quality 
experimental  measurements  and  fully  turbulent 
viscous  flow,  which  would  yield  the  best  conditions 
for  agreement  between  computation  and  experiment. 

Case  4.  Mach  =  3.5,  a  =  8°,  Re/d  =  1.123x10®, 
Figure  5. 

The  laminar  viscous  results  indicate  separation  for 
x/d  =  4.5  and  greater.  Good  agreement  is  shown  on 
the  wind-side  for  all  stations  except  at  x/d  =  2.4, 
where  the  calculations  predict  a  higher  pressure  than 
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the  experiment.  No  indication  of  a  pressure  jump 
resulting  from  an  upstream  shock  is  apparent  as  was 
seen  in  the  lower  Mach  number  flows.  The  turbulent 
viscous  predictions  for  x/d  >7.5  predict  a  lower  first 
suction  peak  than  indicated  by  the  experimental 
data. 

Case  5.  Mach  =  3.5,  a  =  14  Re/d  =  1.123xl0^ 
Figure  6. 

The  laminar  and  turbulent  viscous  results  are  close 
agreement  over  the  full  missile  body.  The  separated 
flow  region  on  the  lee-side  is  characterized  by  a 
relatively  flat  pressure  distribution  compared  to  the 
other  test  cases.  For  x/d  >  5.5,  the  turbulent  results 
indicate  earlier  separation  than  the  laminar  results 
with  the  experiment  located  in  between. 

Case  6.  Mach  =  0.7,  a  =  14  °,  Re/d  =  0.667xl0^ 
Figure  7. 

This  subsonic  test  case  was  added  during  the  second 
year  of  the  study  as  a  special  challenge. 
Computational  modeling  required  moving  the  outer 
boundary  further  away  from  the  body  in  order  to 
lessen  effects  of  the  outer  boundary  interacting  with 
the  flow  on  the  body.  This  reduced  the  density  of  the 
grid  resolution  and  was  cause  for  some  concern.  The 
turbulent  viscous  predictions  are  obviously  in  better 
agreement  with  the  experiment  than  laminar.  The 
lee-side  suction  peak  is  not  predicted  well  until  x/d 
>  7.5  is  reached.  Also,  the  turbulent  results  indicate 
some  jitter  in  the  surface  pressure  profiles, 
particularly  on  the  wind-side  for  x/d  >  3.5. 

In  summary,  the  test  case  results  indicate  that 
turbulent  viscous  modeling  is  required  and  that  the 
flow  fields  are  highly  separated  on  the  lee-side. 

There  is  some  indication  that  a  wind-side  shock 
occurs  resulting  in  a  pressure  jump  on  the  wind-side 
at  downstream  stations  for  test  cases  with  free  stream 
Mach  Number  <  2.0  that  is  not  predicted  by  the 
computations. 

Grid  Resolution 

The  question  of  grid  resolution  was  addressed  in  this 
study  by  making  comparisons  for  several  grid 
resolutions.  Since  the  time-marching  codes  require 
substantial  computer  resources  for  modest  grid 
resolution,  this  effect  has  been  addressed  by 
comparing  results  for  time-marching  codes  using  a 
variety  of  grid  resolutions  and  also  by  comparing  the 
results  from  time-marching  codes  with  much  higher 


resolution  results  obtained  with  an  iterative  PNS 
code.  An  example  comparing  results  for  different 
grid  resolutions  using  the  time-marching  code 
CFL3D  is  shown  in  Figure  8.  Only  slight 
differences  are  seen  for  the  grid  resolutions  shown 
(61  X  81  X  141  and  81  X  111  x  141)(axial,  normal, 
circumferential).  Grid  resolutions  for  the  time¬ 
marching  codes  compared  in  the  study  ranged  from 
61  to  150  (axial),  81  to  lll(normal),  and  81  to 
141  (circumferential) . 

The  iterative  PNS  technique  provides  the  capability 
to  achieve  highly  dense  grid  resolution  for  the 
conditions  of  this  study  with  moderate  computer 
resources.  An  example  of  the  results  achieved  using 
this  technique  for  a  grid  resolution  of  1000  x  72  x 
72  is  shown  in  Figure  9.  The  turbulent  viscous 
result  is  comparable  to  that  shown  in  Figure  8. 

These  results  lend  confidence  that  a  nominal  grid 
resolution  of  80  x  80  x  80  is  adequate  for  the  test 
cases  considered  here.  Figure  9  also  shows  results 
for  laminar  and  inviscid  modeling.  The  inviscid 
result  predicts  a  lee-side  shock  that  is  not  confirmed 
by  the  experimental  measurements. 

Turbulence  Models 

A  variety  of  turbulence  models  have  been  utilized  for 
the  comparisons  with  experimental  data.  The  listing 
includes:  Baldwin-Lomax(BL);  Baldwin-Lomax- 
Degani-Schiff(BLDS);  k-epsilon(ke),  k-omega(kw), 
Spalart-Allmaras(SA);  Smagorinsky(SM);  and 
Baldwin-Barth(BB).  An  example  of  comparisons  for 
Case  3  using  the  NPARC  code  for  four  turbulence 
models  is  shown  in  Figure  10.  The  BB  model 
appears  to  provide  the  best  results  for  the  first 
suction  peak.  The  BL  model,  where  the  turbulent 
length  scale  is  determined  by  a  search  to  the  flow 
field  outer  boundary,  consistently  overpredicts  the 
location  of  the  first  suction  peak.  None  of  the 
models  predict  well  the  separated  flow  region  for 
x/d  =  6.5  and  7.5,  where  there  is  a  distinct  second 
suction  peak.  For  x/d  =  9.5  and  1 1.5,  the  models 
provide  equivalent  results  in  the  lee-side  separated 
flow.  These  results  are  consistent  with  those  for  all 
the  computational  codes  utilized.  The  message 
indicated  here  is  that  the  more  computer-intensive 
two-equation  turbulence  models  do  not  provide  a 
clear  advantage  over  the  less  sophisticated  algebraic 
or  one-equation  turbulence  models. 
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Statistical  Analysis 

The  scope  of  the  computational  study  includes  eight 
computer  codes,  nine  turbulent  viscous  models, 
laminar  viscous  results,  several  grid  resolutions,  and 
six  test  cases.  This  has  resulted  in  more  than  90  test 
case  results  for  evaluation  with  respect  to  each  other 
and  with  respect  to  effects  of  grid  resolution, 
turbulence  model,  and  computational  technique. 

The  traditional  method  for  evaluating  the  previously 
mentioned  effects  on  flow  field  computations  has 
been  for  the  engineer  to  visually  inspect  the  results, 
make  qualitative  comparisons,  and  draw 
conclusions.  This  process  can  be  very  satisfactory 
for  limited  data  sets  to  evaluate.  However,  for  this 
study,  the  number  of  data  sets  has  grown  to  a  large 
number,  thus  making  the  evaluation  and  comparison 
of  the  results  a  formidable  task. 

Thus,  it  was  decided  to  explore  the  possibility  that 
statistical  analysis  techniques  could  assist  in  the 
evaluation  process.  The  technique  used  here  is  to 
consider  the  difference  between  the  computation  and 
experiment.  This  difference  has  been  used  to  obtain 
a  mean  value,  standard  deviation,  and  variance 
between  the  experiment  and  computation.  An 
example  of  the  results  of  this  analysis  technique  for 
Case  3  in  which  the  statistical  quantities  have  been 
summed  for  41  circumferential  positions  at  9  axial 
stations  for  each  computational  result  is  shown  in 
Figure  1 1.  A  visualization  of  the  analysis  is  shown  in 
Figures  1 1  and  12  in  the  form  of  a  box  and  whisker 
plot.  The  summary  shown  consists  of  the  median  of 
the  errors,  the  inner-quartiles  (the  25th  and  75th 
percentiles)  that  determine  the  vertical  dimension  of 
the  box,  and  the  maximum  and  minimum  error  to 
determine  the  whisker  lengths.  These  results 
indicate  that,  although  individual  comparisons  for 
different  codes,  grid  resolutions,  and  viscosity 
models  show  varying  amounts  of  disagreement,  the 
overall  quality  of  the  comparisons  with  experiment  is 
comparable. 

The  above  analysis  technique  has  been  applied  to  all 
data  sets  from  the  study.  The  results  shown  in 
Figures  1 1  and  12  provide  an  overview  of  the 
comparisons  obtedned  for  the  various  codes  and 
permit  some  general  conclusions  to  be  made. 

All  turbulent  viscous  results  for  test  case  3  are  shown 
in  Figure  11.  This  provides  a  good  summary 
comparison  for  the  different  computational 
techniques  and  turbulence  models.  Table  4  provides 


a  listing  of  the  CFD  code  associated  with  the  ID 
number  shown  in  Figure  11.  The  results  indicate 
relative  consistency.  Results  number  1,  4,  and  13 
appear  to  have  the  best  agreement,  whereas  results 
number  2,  8,  10,  12,  and  16  appear  to  have  the  larger 
errors.  Of  these  results,  case  4  and  13  were  obtained 
using  the  BLDS  turbulence  model  and  case  2,  12,  16 
were  obtained  using  the  BL  turbulence  model 
without  modification.  These  results  suggest  that  the 
BL  turbulence  model  has  potential  for  good 
performance;  however,  it  requires  a  skilled 
computational  specialist  to  use  effectively.  Another 
observation  is  that  the  benefit  of  using  more  complex 
two-equation  turbulence  models  does  not  justify  the 
increased  computational  effort  as  opposed  to 
algebraic  or  single  equation  turbulence  models. 

A  sample  statistical  analysis  result  is  shown  in 
Figure  12,  which  shows  results  achieved  using  the 
FDL3DI  code  for  test  cases  1,  2,  3,  4  and  5.  This 
plot  provides  insight  into  the  relative  accuracy  for 
the  predictions  for  the  various  flow  field  conditions. 
One  obvious  observation  is  that  the  results  for  test 
cases  1,  2,  and  3  (Mach  <  3)  indicate  significantly 
greater  errors  than  those  for  test  cases  4  and  5 
(Mach  =3.5). 

This  technique  shows  promise  for  providing  a  useful 
tool  for  use  by  the  engineer  when  comparing 
numerous  computational  and  experimental  results; 
however,  it  does  not  diminish  the  value  of 
examination  of  individual  profile  comparisons. 

Vortical  Flow 

The  experimental  data  included  measurements  of 
pitot  pressures  in  the  outer  flow  fields  using  total 
head  probes.  These  data  provide  the  opportunity  to 
evaluate  the  ability  to  predict  the  size,  strength,  and 
location  of  these  vortices.  These  comparisons  are 
made  by  direct  comparisons  of  pitot  pressures 
calculated  from  the  flow  field  computational  data  to 
the  experimental  measurements.  Additionally,  in  an 
effort  to  obtain  a  more  quantitative  measurement  of 
the  ability  of  the  computational  technique  to  predict 
the  strength  and  location  of  the  vortex,  the  location 
of  the  center  of  the  vortex  and  the  pitot  pressure 
along  a  line  extending  from  the  vortex  center  to  the 
pitch-plane  axis  were  obtained  from  the  experiment. 

Examples  of  comparisons  between  the  computations 
and  experiment  for  the  pitot  pressure  surveys  are  shown 
in  Figures  13, 14, 15,  and  16  for  M=  1.8,  2.5, 3.5,  and 
0.7,  respectively  for  x/d  =  1 1 .5.  These  figures  show 
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examples  of  pitot  contours  and  the  pitot  pressure 
extrapolation  along  a  line  perpendicular  to  the  pitch 
plane  passing  through  the  vortex  core  as  defined  by  the 
experimental  measurement.  Comparisons  are  shown 
for  laminar  and  turbulent  viscous  modeling  obtained 
using  the  NPARC  code  and  the  BLDS  turbulence 
model.  The  results  indicate  that  the  laminar  viscous 
predictions  consistently  underpredict  the  strength  and 
size  of  the  vortex  core.  Also,  the  turbulent  modeling 
provides  encouraging  accuracy  for  the  axial  position 
shown. 

The  NPARC  implementation  of  the  BLDS  turbulence 
model  provides  an  option  to  limit  the  search  range  in 
the  radial  direction  for  determining  the  turbulent  length 
scale.  An  example  of  the  effect  of  this  limiter  on  the 
prediction  of  the  outer  vortical  flow  is  shown  in  Figure 
17.  The  case  shown  is  for  M=2.5,  (X^14°.  Results  for  a 
search  limit  of  15  and  25  points  are  shown  (ke  =  15  and 
ke=25).  The  extent  of  the  region  of  this  search  in  the 
flowfield  is  indicated  by  the  circles  drawn  about  the 
model  cross  section.  A  detailed  comparison,  not  shown, 
in  the  study  was  carried  out  to  examine  this  effect  on  the 
predictions.  It  was  determined  that  for  the  test  cases 
considered  here,  a  value  of  15  provided  optimal  results. 

Aerodynamic  Forces  and  Moments 

The  experimental  data  available  included  strain  gage 
balance  measurements  for  normal  force,  pitching 
moment,  and  axial  force.  A  partial  tabulation  of 
results  obtained  for  normal  force  and  pitching 
moment  are  provided  in  Table  5.  These  results 
indicate  that  most  codes  provided  comparable 
predictions  that  are  within  ±3%  agreement  with 
experiment  for  normal  force  and  ±  5%  for  pitching 
moment.  This  degree  of  accuracy  is  within 
acceptable  requirements  for  initial  design  and 
optimization  analysis  for  most  projectile  and  missile 
systems. 

Closing  Remarks 

1)  Consistent  results  were  achieved  for  the  various 
computational  codes  for  like  conditions  of  grid 
density  and  turbulence  model. 

2)  For  the  conditions  of  this  study,  no  best  turbulence 
model  could  be  identified  in  the  comparisons  to 
experiment.  There  was  no  consistent  benefit  obtained 
from  utilization  of  a  (k-epsilon  or  k-omega) 
turbulence  model  over  the  algebraic  models  even 


though  substantial  flow  separation  was  present  in  the 
flow  fields  studied. 

3)  Although  flow  field  details  such  as  the  definition 
of  the  separated  vortex  and  the  surface  pressure  in 
regions  of  separated  flow  are  not  predicted  as  well  as 
desired,  the  aerodynamic  force  and  moment 
predictions  are  accurate  within  ±  5%  for  these  test 
cases.  This  level  of  accuracy  is  usually  adequate  for 
weapon  design  requirements.  Thus,  CFD  is  a  viable 
predictive  tool  for  the  class  of  body  configurations 
considered  in  this  study  and  should  be  of  increasing 
importance  as  computer  technology  continues  to 
advance  thus  permitting  the  predicted  results  to  be 
obtained  in  a  timely  fashion. 
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Nomenclature 

d  model  cylinder  diameter,  3.7-inches 

Cm  pitching  moment  coefficient 

Cn  Normal  force  coefficient 

M  Mach  Number 

X  axial  distance  from  nose  of  missile 
y  vertical  distance  from  axis 

z  horizontal  distance  from  axis 

a  angle  of  attack,  degrees 

phi  circumferential  angle  about  axis 

Turbulence  Model  Identification 

BB  Baldwin-Barth 

BBGR  Baldwin-Barth-Goldstein- 
Ramakrishnan 
BL  Baldwin-Lomax 

BLDS  Baldwin-Lomax  with  Degani-Schiff 
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ke  k-epsilon 

kw  k-omega,  Wilcox 

SA  Spalart-Allmaras 

SM  Smagorinski 

kwm  k-omega,  Menter 
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Table  2.  Wind  Tunnel  Test  Conditions 


Case  Number 

Mach  Number 

Reynolds  No,  xlO^/ft 

1 

1.45 

14 

2.0 

2 

1.8 

14 

2.0 

3 

2.5 

14 

4.0 

4 

3.5 

8 

4.0 

5 

3.5 

14 

2.0 

6 

0.7 

14 

2.0 

6 


Table  3.  Experimental  Data  Available 


Case  Number 

Surface  Pressure 

Pitot  Survey,  X/D 

Force  Measurements 

1 

32x41 

8.5, 11.5 

no 

2 

32x41 

5.5, 8.5, 11.5 

no 

3 

32x41 

5.5, 11.5 

yes 

4 

32x41 

5.5 

yes 

5 

32x41 

5.5, 11.5 

yes 

6 

32x41 

11.5 

no 

Table  4.  Code  ID  for  Figure  11  and  Table  5. 


ID  Number 

Code,  Turbulence  Model 

1 

AHPCRCFE,SM 

2 

UNPNS,BL 

3 

UNIS/USA,  BBGR 

4 

OVERFLOW,  BEDS 

5 

FDL3DI,  ke 

6 

COBALT,  SA 

7 

NPARC,BLDS 

8 

CFL3D,  BB 

9 

CFL3D,  BL 

10 

CFL3D,kw 

11 

CFL3D,SA 

12 

NPARC,  BB 

13 

NPARC,BL 

14 

NPARC,  ke 

15 

NPARC,  kw 

16 

F3D,  BL 

17 

USA,BB 

Table  5.  Normal  Force  and  Pitching  Moment  -  Prediction  and  Experiment 


Figurel.  Schlieren  photograph  of  model  mounted  in  wind  tunnel,  Mach  =3.5,  a  =  8®. 
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Figure  3.  Test  Case  Overview,  Case  2,  Mach  =  1.8,  a  =  14°,  Re/D  =  0.667x10® 
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Figure  5.  Test  Case  Overview,  Case  4,  Mach  =  3.5,  a  =  8°,  Re/D  =  1.123x10^ 
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Figure  7.  Test  Case  Overview,  Case  6,  Mach  =  0.7,  a  =  14°,  Re/D  =  0.667x10 
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Figure  8.  Grid  resolution,  CFL3D  results 
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Figure  9.  Grid  resolution,  iterative  PNS,  Case  3,  Mach  =  2.5,  a  =  14°,  X/D  =  1 1.5 
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ire  13.  Pitot  pressure  contours  for  experiment  cc 
iputations.  Vortex  core  pitot  pressure  along  Z/D 
X/D=11.5. 
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Figure  14.  Pitot  pressure  contours  for  experiment  compared  to  laminar  and  turbulent  viscous 
computations.  Vortex  core  pitot  pressure  along  Z/D  line  to  flow  field  centerline. 

Case  3,  M=2.5, 0^14°,  X/D  =1 1.5. 


Figure  1 6.  Pilot  pressure  contours  for  experiment  compared  to  laminar  and  turbulent  viscous 
computations.  Vortex  core  pitot  pressure  along  Z/D  line  to  flow  field  centerline. 

Case  6,  M=0.7,  oc^l4°,  X/D  =1 1.5. 


APPENDIX  A.2 

Technical  Paper,  “Structured  Grid  Based  Solution-Adaptive 
Technique  for  Complex  Separated  Flows”,  to  appear. 
Journal  of  Applied  Mathematics  and  Computation 
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A  Structured  Grid  Based  Solution-Adaptive  Technique  for  Complex  Sepa- 

■  ■  •  rated  Flows  • 

by 

Hugh  Thornburg,  Bharat  K.  Soni  and  Boyalakuntla  Kishore 
NSF  Engineering  Research  Center  for 
Computational  Field  Simulation 
Mississippi  State  University 
Mississippi  State,  MS  39762 

ABSTRACT 

A  structured  grid  based  technique  is  presented  to  enhance  the  predictive  capability  of  widely  used 
CFD  codes  through  the  use  of  solution  adaptive  gridding.  The  procedure  redistributes  mesh  points 
based  upon  an  existing  flowfield.  A  newly  developed  weight  function  employing  Boolean  sums  is 
utilized  to  represent  the  local  truncation  error.  This  weight  function  is  then  used  to  construct  forcing 
functions  associated  with  an  elliptic  system!  A  NURBS  representation  is  employed  to  define  block 
surfaces  for  boundary  point  redistribution.  The  technique  has  been  applied  to  a  flowfield  about  a 
configuration  of  practical  interest.  This  flowfield  involves  supersonic  fir eestream  conditions  at  angle 
of  attack,  and  exhibits  large  scale  separated  vortical  flow,  vortex-vortex  and  vortex-surface  interac¬ 
tions,  separated  shear  layers  and  multiple  shocks  of  different  intensity.  The  results  demonstrate  the 
capability  of  the  developed  weight  function  to  detect  shocks  of  differing  strengths,  primary  and  sec¬ 
ondary  vortices,  and  shear  layers  adequately. 

INTRODUCTION 

Rapid  access  to  highly  accurate  data  about  complex  configurations  is  needed  for  multi-disciplinary 
optimization  and  design.  In  order  to  efficiently  meet  these  requirements  a  closer  coupling  between 
the  analysis  algorithms  and  the  discretization  process  is  needed.  In  some  cases,  such  as  free  surface, 
temporally  varying  geometries,  and  fluid  structure  interaction,  the  need  is  unavoidable.  In  other 
cases  the  need  is  to  rapidly  generate  and  modify  high  quality  grids.  Techniques  such  as  unstructured 
and/or  solution-adaptive  methods  can  be  used  to  speed  the  grid  generation  process  and  to  automati¬ 
cally  cluster  mesh  points  in  regions  of  interest.  Global  features  of  the  flow  can  be  significantly  af¬ 
fected  by  isolated  regions  of  inadequately  resolved  flow.  These  regions  may  not  exhibit  high  gradi¬ 
ents  and  can  be  difficult  to  detect.  Thus  excessive  resolution  in  certain  regions  does  not  necessarily 
increase  the  accuracy  of  the  overall  solution. 

Several  approaches  have  been  employed  for  both  structured  and  unstructured  grid  adaption.  The 
most  widely  used  involve  grid  point  redistribution,  local  grid  point  emichment/derefinement  or  lo¬ 
cal  modification  of  the  actual  flow  solver.  However,  the  success  of  any  one  of  these  methods  ulti¬ 
mately  depends  on  the  feature  detection  algorithm  used  to  determine  solution  domain  regions  which 
require  a  fine  mesh  for  their  accurate  representation.  Typically,  weight  functions  are  constructed  to 
mimic  the  local  truncation  error  and  may  require  substantial  user  input.  Most  problems  of  engineer¬ 
ing  interest  involve  multi-block  grids  and  widely  disparate  length  scales.  Hence,  it  is  desirable  that 
the  adaptive  grid  feature  detection  algorithm  be  develojjed  to  recognnse  flow  structures  of  different 
type  as  well  as  differing  intensity,  and  adequately  address  scaling  and  normalization  across  blocks. 
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These  weight  functions  can  then  be  used  to  constnict  blending  functions  for  algebraic  redistribution, 
interpolation  functions  for  unstructured  grid  generation,  forcing  functions  to.  attract/tepel  points  in 
an  elliptic  system,  or  to  trigger  local  refinement,  based  upon  application  of  an  equidistribution  prin¬ 
ciple.  The  popularity  of  solution-adaptive  techniques  is  growing  in  tandem  with  unstructured  meth¬ 
ods.  The  difficultly  of  precisely  controlling  mesh  densities  and  orientations  with  current  unstruc¬ 
tured  grid  generation  systems  has  driven  the  use  of  solution-adaptive  meshing  Use  of  derivatives 
of  density  or  pressure  are  widely  used  for  construction  of  such  weight  functions,  and  have  been  prov¬ 
en  very  successful  for  inviscid  flows  with  shocks[2,7,ll].  However,  less  success  has  been  realized 
for  flowfields  with  viscous  layers,  vortices  or  shocks  of  disparate  strength.  It  is  difficult  to  maintain 
the  appropriate  mesh  point  spacing  in  the  various  regions  which  require  a  fine  spacing  for  adequate 
resolution.  Mesh  points  often  migrate  from  important  regions  due  to  refinement  of  dominant  fea¬ 
tures.  An  example  of  this  is  the  well  know  tendency  of  adaptive  methods  to  increase  the  resolution 
of  shocks  in  the  flowfield  around  airfoils,  but  in  the  incorrect  location  due  to  inadequate  resolution 
of  the  stagnation  region.  This  problem  has  been  the  motivation  for  this  research. 

The  weight  functions  developed  utilize  scaled  derivatives  and  normahzing  procedures  to  minimize 
or  eliminate  the  need  for  user  input.  The  most  attractive  feature  of  this  work  is  the  ability  to  detect 
flow  features  of  varying  intensity  and  the  lack  of  user  defined  inputs  for  the  selection  of  the  weight 
function. 

In  this  research  a  NURBS  representation  is  employed  to  define  block  surfaces  for  boundary  point 
redistribution.  The  features  described  have  been  implemented  into  Adapt2D/3D.  An  adaptive  grid 
system  capable  of  automatically  resolving  complex  flows  with  shock  waves,  expansion  waves,  shear 
iayers  and  complex  vortex-vortex  and  vortex-surface  interactions.  An  adaptive  grid  approach 
seems  well  suited  for  such  problems  in  which  the  spatial  distribution  of  length  scales  is  not  known 
a  priori.  •  . 

APPROACH  TO  ADAPTION 

The  elliptic  generation  system; 

=0  (1) 

1=1  y=i  k=i 

where  r  :  Position  vector, 

g'j  :  Contravariant  metric  tensor 
:  Curvilinear  coordinate,  and 
Pjc  ;  Control  function. 

is  widely  used  for  grid  generation  [  1  ] .  Control  of  the  distribution  and  characteristics  of  a  grid  system 
can  be  achieved  by  varying  the  values  of  the  control  functions  Pk  in  Equation  1  [1].  The  application 
of  the  one  dimensional  form  of  Equation  1  combined  with  equidistribution  of  the  weight  function 
results  in  the  definition  of  a  set  of  control  functions  for  three  dimensions. 


Pi  = 


O'  =  1,2.3) 


These  control  functions  were  generalized  by  Eiseman  [2]  as: 
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In  order  to  conserve  some  of  the  geometric^  Characteristics  of  the  original  grid  the  definition  of  the 
control  functions  is  extended  as:' 
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where  P  initial  geomciry  •  Control  fimction  based  on  initial 

geometry 

Pwt  :  Control  function  based  on 

current  solution 

Cj  :  Constant  weight  factor. 

These  control  functions  are  evaluated  based  on  the  current  grid  at  the  adaption  step.  This  can  be 
formulated  as; 

pM  =  />(«-!)  +  (/•  =  1.2,3)  (5) 

where 

Pf  =  Pf^  +  (/  =  1.2.3)  (6) 

A  flow  solution  is  first  obtained  with  an  initial  grid.  Then  the  control  functions  Pi  are  evaluated  in 
accordance  with  Equations  2  and  5,  which  is  based  on  a  combination  of  the  geometry  of  the  current 
grid  and  the  weight  functions  associated  with  the  cuirent  flow  solution!  11]. 

Evaluation  of  the  forcing  ftmctions  corresponding  to  the  grid  input  into  the  adaptation  program  has 
proven  to  be  troublesome.  Direct  solution  of  Equation  1  for  the  forcing  functions  using  the  input 
grid  coordinates  via  Cramer’s  rule  or  IMSL  libraries  was  not  successfuL  For  some  grids  with  very 
high  aspect  ratio  cells  and/or  very  rapid  changes  in  cell  size,  the  forcing  functions  became  very  large. 
Tlie  use  of  any  differerioing  scheme  other  than  the  one  used  to  evaluate  the  metrics,  sudh'  as  the  hybrid 
upwind  scheme[8],  would  result  in  very  large  mesh  point  movements.  An  alternative  technique  for 
evaluating  the  forcing  functions  based  on  derivatives  of  the  metrics  was  implemented[3]. 
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This  technique  has  proven  to  be  somewhat  more  robust,  but  research  efforts  are  continuing  in  this 
area. 


SOLUTION  PROCEDURE 

The  software  developed  inputs  two  PLOT3D  format  files  [4],  one  for  the  grid  and  one  for  the  flow 
solution.  These  files  contain  the  number  of  blocks ,  the  grid  size,  the  grid  coordinates  and  the  solution 
vector.  Output  consists  of  an  NPARC[5]  restart  file,  as  well  as  two  PLOT3D  files  of  the  adapted 
grid  and  the  flow  solution  interpolated  onto  the  new  grid.  A  multi-block  treatment  c^ability  is  in¬ 
cluded.  The  adaptive  grid  is  constructed  in  three  steps.  The  first  step  is  to  generate  the  weight  func¬ 
tions,  which  due  to  their  critical  importance  will  be  discussed  in  detail  in  a  separate  section.  The 
second  step  is  to  gen^ate  the  actual  adapted  grid  by  equidistribution  of  the  aforementioned  weight 
function.  In  the  current  work  this  is  accomplished  by  the  numerical  solution  of  Equation  1.  A 
coupled  three-dimensiopal  strongly  implicit  procedure  (CSIP)  as  desaibed  by  Ghia  et  al.  [6]  has 
been  implemented  for  the  solution  of  the  discretized  equations.  Upwind  differencing,  with  biasing 
based  on  the  sign  of  the  forcing  functions,  as  well  as  central  differencing  has  been  implemented  and 
studied  for  the  first  derivative  terms.  The  first  order  upwind  differencing  increases  the  stability  of 
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the  procedure,  but  at- the  expense  of  smearing  the  grid  clustering.  Hence,  a  hybrid  upwind/central 
■.differencing  scheme  has  been  implemented  to  lessen  this  smearing  while  maintaining'the  smooth¬ 
ness  and  stability  of  the  upwind  procedure.  Central  differencing  has  been  employed  for  all  second 
and  mixed  derivative  terms.  All  non-linear  terms  were  treated  by  quasi-Unearization.  Since  Equa¬ 
tion  1  is  solved  iteratively,  the  forcing  functions  must  be  evaluated  for  each  new  successive  grid  point 
location.  The  forcing  functions  are  obtained  by  inteipolation  from  the  original  grid.  The  interpola¬ 
tion  procedure  employed  was  that  used  for  the  non-matching  block  to  block  interface  capability  of 
the  NAPRC  [5]  code.  This  scheme  is  based  on  trilinear  interpolation.  Upon  convergence  of  this 
process  the  flow  solution  is  interpolated  onto  the  adapted  grid  using  the  same  interpolation  subrou¬ 
tine.  Calculations  can  then  be  continued  from  the  new  restart  file.  This  procedure  can  then  be  re¬ 
peated  until  an  acceptable  solution  is  obtained.  Experience  indicates  that  coupling  this  procedure 
with  a  code  capable  of  treating  time  accurate  grid  movement  would  ease  this  process  and  lessen  the 
CPU  requirements. 

Currently,  grid  adaptation  is  performed  in  an  uncoupled  manner.  For  simple  problems  this  is  ade¬ 
quate  as  only  marginal  changes  are  observed  after  two  to  three  cycles  of  adaptation.  However,  for 
complex  flow  fields  it  appears  necessary  to  closely  couple  the  flow  solver  and  the  adaptation  tech¬ 
nique  as  significant  changes  in  both  grid  and  flow  solution  are  visible  after  even  six  cycles  of  un¬ 
coupled  adaption.  Also,  the  more  complex  flows  require  the  use  of  a  hybrid  differencing  scheme 
for  the  solution  of  the  grid  equations.  For  flows  containing  strong  shocks  and  shear  layers  along  with 
weaker  structures,  central  differencing  was  unstable  and  first  order  upwind  differencing  of  the  grid 
equations  smeared  the  weaker  structures.  Therefore  to  obtain  full  benefit  of  the  adaptive  grid 
blended  central/upwind  has  been  implemented  for  the  grid  equations. 

WEIGHT  FUNCTIONS 

Application  of  the  equidistribution  law  results  in  grid  spacing  inversely  proportional  to  the  weight 
function,  and  hence,  the  weight  function  determines  the  grid  point  distribution.  Ideally,  the  weight 
would  be  the  local  truncation  error  ensuring  a  uniform  distribution  of  error.  However,  evaluation 
of  the  higher-order  derivatives  contained  in  the  truncation  error  from  the  available  discrete  data  is 
progressively  less  accurate  as  the  order  increases  and  is  subject  to  noise.  Determination  of  this  func¬ 
tion  is  one  of  the  most  challenging  areas  of  adaptive  grid  generation.  Lower-order  derivatives  must 
be  non-zero  in  regions  of  wide  variation  of  higher-order  derivatives,  and  are  proportional  to  the  rate 
of  variation.  Therefore,  lower-order  derivatives  are  often  used  to  construct  a  weight  function  as  a 
proxy  for  the  truncation  error.  Construction  of  these  weight  functions  often  requires  the  user  to  spec¬ 
ify  which  derivatives  and  in  what  proportion  they  are  to  be  used.  This  can  be  a  time  consuming  pro¬ 
cess.  Also,  due  to  the  disparate  strength  of  flow  features,  important  features  can  be  lost  in  the  noise 
of  dominant  features.  The  weight  functions  developed  by  Soni  and  Yang  [7]  and  Thornburg  and  S  oni 
[8]  are  examples  of  such  efforts.  The  weight  function  of  Thornburg  and  Soni  [8]  has  the  attractive 
feature  of  requiring  no  user  specified  inputs.  Relative  derivatives  are  used  to  detect  features  of  vary¬ 
ing  intensity,  so  that  weaker,  but  important  structures  such  as  vortices  are  accurately  reflected  in  the 
weight  function,  fri  addition,  each  conservative  flow  variable  is  scaled  independently.  One-sided 
differences  are  used  at  boundaries ,  and  no-slip  boundaries  require  special  treatment  since  the  veloc¬ 
ity  is  zero.  This  case  is  handled  in  the  same  manner  as  zero  velocity  regions  in  the  field.  A  small 
value,  epsilon  in  equation  8,  is  added  to  all  normalizing  quantities.  In  the  present  work  this  weight 
function  has  been  modified  using  the  Boolean  sum  construction  method  of  Soni  [7].  Also,  several 
enhancements  of  an  implementation  nature  have  been  employed.  For  example  epsilon  has  been 
placed  outside  the  absolute  value  operator.  This  eliminated  the  possibility  of  spurious  gradients  in 
the  weight  function  in  regions  where  epsilon  was  nearly  equal  and  opposite  in  sign  to  the  local 
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normalizing  flow  variable.  Also,  the  normalizing  derivatives  have  been  set  to  an  initial  or  minimum 
value  of  ten  percent  of  the'freestream  quantities.  This  alleviates  problems  encountered  in  flows 
without  significant  features  to  trigger  adaption  in  one  or  more  coordinate  directiohs.  Otherwise  a 
few  percent  variation  would  be  normalized  to  the  same  level  as  a  shock  or  other  strong  feature.  The 
current  weight  function  is  as  follows: 

w  -  ffi 

max(W\W\VP)  max(W*,W\W^) 
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The  symbol  ®  represents  the  Boolean  sum.  Note  that  the  directional  weight  functions  are  scaled 
using  a  common  maximum  in  order  to  maintain  the  relative  strength.  For  the  results  shown,  the 
Wdight  fuiifbtioh  used  is  the  sum  of  the  weight  function  defined  by  Eq.  8  and  the  one  defined  in  [8]. 
The  weight  function  defined  in  [8]  is  Eq.  8  evaluated  without  using  the  Boolean  sum. 

Flow  past  a  wedge  at  Mach = 2.0  is  used  to  illustrate  the  enhanced  detection  capabilities  of  this  newly 
developed  weight  function.  Figure  1  presents  weight  functions  evaluated  using  the  previous  proce¬ 
dure,  lower  half  plane,  as  well  as  the  current  procedure,  upper  half  plane. 


ADAPTED  GRID  Wim^TTHOCT (BOOLEAN) 


WEIGHT  WrrmTTHOUT  BOOLEAN 


Figure  1.  Comparison  of  Weight  Functions. 
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Figure  2.  Comparison  of  Adapted  Grids. 


5 


It  can  be  observed  that  both  weight  functions  clearly  detected  the  primary  shock.  It  can  also  be  seen 
that  the  expansion  fan,  boundary  layer,  arid  the  reflect^  shocks  are  much  more  clearly  represented 
in  the  current  weight  function.  Adapted  grids  using  both  weight  function  formulations  are  presented 
in  Fig.  2.  The  high  gradient  regions  of  the  expansion  region  are  only  reflected  in  the  adapted  grid 
using  the  new  weight  function.  The  reflected  shock  is  also  much  sharper.  Figure  3  compares  the 
solution  obtained  using  the  current  adaption  procedure  with  that  obtained  using  the  original  grid. 
The  enhanced  resolution  is  clearly  evident. 
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BOUNDARY  POINT  REDISTRIBUTION 

Accurate  representation  of  the  flowfield  in  the  vicinity  of  boimdaries  is  critical  for  an  acceptable 
overall  solution.  Physics  occurring  near  boundaries  often  drive  the  flow  physics  occurring  in  other 
regions  of  the  flowfield.  This  is  especially  true  for  noslip  surfaces.  Hence,  the  quality  and  distribu¬ 
tion  of  the  grid  in  this  region  is  of  critical  importance.  For  the  implementation  of  many  turbulence 
models  orthogonality  is  also  required.  When  using  an  adaptive  procedure  based  on  a  redistribution 
of  mesh  points,  such  as  in  this  work,  the  interior  points  move  as  the  grid  is  adapted.  This  leads  to 
distorted  cells  if  the  boundary  points  are  not  redistributed  as  the  grid  is  adapted.  Both  grid  quality 
and  geometric  fidelity  must  be  maintained  during  redistribution.  In  the  current  work  aU  surfaces  of 
individual  blocks  are  treated  in  the  same  manner,  whether  block  interfaces  or  flow  boundary  condi¬ 
tions.  Two  steps  must  be  performed.  First  the  geometry  is  defined,  and  secondly  the  surface  mesh 
is  regenerated  using  a  given  distribution  mesh.  The  geometry  is  defined  as  a  NURB  surface  from 
the  current  surface  mesh  to  be  redistributed.  This  procedure  is  based  on  the  algorithms  developed 
by  Yu  and  Soni  for  Genie++[9].  This  definition  is  then  used  to  generate  the  surface  using  a  user 
specified  distribution  mesh.  The  entire  surface  or  a  subregion  can  be  redistributed.  Subregions  can 
be  used  to  fix  points,  such  as  sharp  comers  or  a  transition  point  between  boundary  condition  type. 
For  example  block  boimdary  to  noslip  surface.  The  distribution  mesh  is  a  [0,1]  square  evaluated 
from  a  specified  surface  based  upon  arc  length.  For  solid  surfaces  the  distribution  mesh  is  based 
on  the  nearest  interior  coprdinate  surfaces.  The  spacing  between  surfaces  is  small  and  the  surfaces 
are  of  a  similar  geometric  shape.  Therefore  the  normal  coordinate  is  nearly  orthogonal  to  the  surface. 
Block  interfaces  are  treated  by  redistributing  the  current  block  surface  based  on  its  corresponding 
surface  in  the  neighboring  block. 
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Figure  4.  Adapted  grid  after  two  cycles. 


Figures.  Adapted  grid  after  two  cycles. 


Figure  5  presents  the  grid  constructed  using  the  previous  weight  function  and  the  same  flow  condi¬ 
tions  and  number  of  adaptation  cycles.  Figures  5  and  6  present  streamwise  cuts  of  the  two  grids 
shown  in  Figs  4  and  5  at  X/D  =  5.5  and  7.5  respectively.  The  left  had  side  of  Figures  6  and  7  corre¬ 
spond  to  Figure  5  and  the  right  hand  side  corresponds  to  Figure  5. 


As  can  be  seen  from  Figures  4-7  the  feeding  sheet,  and  both  primary  and  secondary  vortices,  as  well 
as  the  strong  shock  at  the  nose  are  clearly  visible.  Comparison  of  Figures  4  and  5  show  the  sharper 
resolution  provided  by  the  new  formulation  of  the  weight  function.  The  shock  itself  is  more  sharply 
represented,  but  the  greatest  improvement  is  in  representing  the  vortex  and  the  shear  layer.  The  pre¬ 
vious  formulation  simply  clustered  mesh  points  in  the  vicinity  of  the  vortex.  This  does  improve  reso¬ 
lution.  However,  the  new  formulation  clearly  shows  the  circular  shape  of  the  vortex,  as  well  as  the 


high  gradient  regions  on  the  edges  of  the  structure.  Hence,  the  new  method  of  grid  adaptation  more 
closely' reproduces  the  flow  physics.  Also,  stronger  and  sharper  concentration  of  mesh  points  is 
observed  with  the  new  formulation.  This  is  most  apparent  in  the  far  field  region  and  the  greater  detail 
of  the  flowfield  visible  in  the  grid.  Both  the  primary  and  secondary  separation  points  are  well  repre¬ 
sented.  On  the  windward  side  of  the  missile  it  can  be  seen  that  the  boundary  layer  clustering  is  great¬ 
ly  increased  by  the  adaptation  procedure.  This  is  desirable  in  this  area  as  the  boundary  layer  is  the 
only  feature  of  interest.  Examination  of  Fig  4  reveals  that  toward  the  aft  end  of  the  missile,  the  off 
surface  structure  is  not  well  defined.  This  is  because  after  two  cycles  of  adaptation  the  structure  is 
not  well  enough  defined  by  the  flow  solver  and  associated  grid  for  the  adaptation  procedure  to  sharp¬ 
ly  define  the  structure.  Only  structures  that  have  been  at  least  partially  resolved  by  the  flow  solver 
can  be  detected  by  the  weight  function.  Hence  the  quality  of  the  weight  function  is  dependent  upon 
the  quality  of  the  solution.  It  should  be  noted  that  resolution  of  the  flowfield  improves  with  each 
cycle  of  adaptation.  The  adaption  procedure  and  the  flow  solver  should  be  coupled  so  that  the 
adapted  grid  can  reflect  all  the  features  that  are  detected  as  the  solution  progresses  and  improves  due 
to  adaption. 

Figure  8  presents  the  flow  solution  obtained  using  the  NPARC  [5]  flow  solver,  and  the  KE  turbu¬ 
lence  model  option.  Figure  9  presents  the  associated  weight  function. 


Figure  8.  Normalized  Stagnation  Pressure.  Figure  9.  Weight  Function. 


CONCLUSIONS 

The  results  shown  demonstrate  the  capability  of  the  developed  weight  function  to  detect  shocks  of 
differing  strengths,  primary  and  secondary  vortices,  and  shear  layers  adequately.  For  the  test  case 
presented  the  increased  resolution  of  the  shock  and  the  expansion  region  resulted  in  the  structure 
increasing  in  size  so  that  interaction  with  the  farfield  occurred.  Thus,  the  farfield  boundaries  should 
be  placed  at  a  larger  distance  from  the  body.  No  user  input  is  required.  It  is  imperative  that  the 
adaptation  process  be  coupled  with  the  flow  solver  as  experience  indicates  that  the  complex  flow- 
fields  may  require  many  adaptive  cycles.  This  is  particularly  critical  for  flow  fields  that  are  not  well 
represented  by  the  initial  solution.  The  adapted  grids  allow  the  use  of  larger  time  steps  to  increase 
the  convergence  rate.  However,  the  single  greatest  benefit  resulting  from  the  adapted  grid  is  the  low¬ 
ering  of  the  artificial  viscosity  required  for  stability.  The  adapted  grid  aligns  and  clusters  near  shocks 
and  shear  layers.  For  the  test  problem  a  dramatic  decrease  in  the  artificial  viscosity  coefficients  is 
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possible  when  running  on  the  adapted  grid.  This  results  in  off  surface  structures  which  are  much 
sharper  and  more  closely  resemble  the  flow  visualization.  The  boundary  layer  on  the  windward  side 
also  became  thinner  and  the  normal  grid  spacing  was  decreased  by  the  adaptive  procedure. 

Currently  multi-block  problems  are  being  simulated  to  evaluate  this  capability.  Attention  is  being 
placed  upon  across  block  scaling  for  the  weight  functions.  A  global  maximum  across  all  blocks 
seems  to  work  well  for  the  normalization  of  weight  function  components.  However,  this  issue  is  be¬ 
ing  monitored  as  more  complex  problems  are  attempted.  A  further  problem  arises  due  to  equidis- 
tribution  via  forcing  functions.  This  does  not  appear  to  be  a  problem  where  a  block  face  is  connected 
to  one  and  only  one  face  of  another  block.  Problems  have  been  encountered  for  a  multi-block  launch 
vehicle  computation  where  a  block  face  consisted  of  a  block  to  block  connection  and  a  solid  surface 
for  the  base  region.  A  smooth,  continuous  set  of  forcing  functions  across  this  block  boundary  does 
not  yield  a  smooth  grid.  It  is  believed  that  this  is  due  to  the  elliptic  character  of  the  grid  equations. 
The  simple  fix  is  to  introduce  artificial  block  boundaries  to  force  one  to  one  correspondence.  Work 
is  continuing  on  these  issues  as  well  as  coupling  with  a  flow  solver,  unsteady  flows,  and  more  com¬ 
plex  multi-block  problems.  These  results  will  be  published  at  a  later  date. 
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